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Abstract 

The shear viscosity of a two-dimensional (2D) liquid was calculated using equilibrium molecular 
dynamics simulations with a Yukawa potential. The shear viscosity has a minimum, at a Coulomb 
coupling parameter V of about 17, arising from the temperature dependence of the kinetic and 
potential contributions. Previous calculations of 2D viscosity were less extensive, and for a different 
potential. The stress autocorrelation function was found to decay rapidly, contrary to some earlier 
work. These results are useful for 2D condensed matter systems and are compared to a recent 
dusty plasma experiment. 

PACS numbers: 52.27.Lw, 82.70.Dd, 52.27.Gr 
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Two-dimensional systems in crystalline or liquid states 1] are of interest in various fields 
of physics. Monolayer particle suspensions can be formed in colloidal suspensions 3] and 
dusty plasmas jjj. Electrons on the surface of liquid helium form a 2D Wigner crystal 0. 
Ions in a Penning trap can be confined as a single layer of a one-component plasma (OCP) |5| . 
Magnetic flux lines in 2D high-temperature superconductors form patterns of hexagonally 
correlated vortices 61. At an atomic scale, gas atoms adsorb on the surface of substrates 

Q 

such as graphite f7jj. Here we will be concerned with liquids, including liquids near freezing, 
composed of molecules or particles that interact with a Yukawa pair potential. 



rese include colloids, mono- 
in biological and chemical 



The Yukawa pair potential is widely used in several fields. T 
layer strongly-coupled dusty plasmas, and some polyelectrolytes 
systems. The Yukawa potential energy for two particles of charge Q separated by a distance 
r is U(r) = Q 2 (47reor) _1 exp(— r/Xn), where Xd is a screening length. This potential changes 
gradually from a long-range r~ l Coulomb repulsion to a hard-sphere- like repulsion as the 
screening parameter k = a/Xo is increased. Here, a = (mr) -1 / 2 is the 2D Wigner-Seitz 
radius [9] and n is the areal number density of particles. 



The literature has only a few reports o 



puting the shear viscosity of 2D liquids 



10, 



molecular dynamics (MD) simulations for corn- 



Ill , and none of those are for a Yukawa potential. 



Here we present such a simulation, yielding results for the shear viscosity and the shear stress 
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autocorrelation function [SACFY. We will compare results to a recent 2D experiment 
and to MD simulations |l3j, llJ, lla, llfij for the shear viscosity of 3D liquids. 

Our first motivation arises from the need to model the recent dusty plasma experiment 
of Nosenko and Goree jl^. The experiment was performed with a monolayer of polymer 
microspheres suspended in a plasma. The microspheres interacted with a Yukawa pair po- 



tential 17[. In an undisturbed state, particles arranged themselves in a 2D triangular lattice, 
which was then melted by an externally-applied velocity shear due to two counterpropagat- 
ing laser beams applied in situ. In this way, a 2D liquid was produced that had a shear flow. 
The experimenters measured rj and found its variation with temperature. 

Our second motivation is to compare rj for 2D and 3D liquids, both with a Yukawa 
potential. Because shear viscosity has the different units of kg m _1 s _1 and kg s" 1 in 3D 
and 2D, respectively, we divide by the volume and areal mass density respectively, yielding 
the kinematic viscosity v. This quantity has the same units of m 2 s _1 for both 3D and 2D, 
thereby allowing a comparison of results for 3D and 2D. 
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Our simulation uses an equilibrium method to calculate r\. Under equilibrium conditions, 
momentum transport arises from random thermal fluctuations in a homogeneous sample, 
and there is no macroscopic shear flow. In contrast, we note that in a shear flow the system 
is not in equilibrium, and the viscosity often depends on the applied shear rate. 

In an equilibrium method, shear viscosity can be calculated using the Green-Kubo rela- 
tion [18]. Green-Kubo formulae in general yield a macroscopic phenomenological transport 
coefficient, such as the viscosity or diffusion coefficient, which is written as a time integral of 
a microscopic time-correlation function. The Green-Kubo approach is based on a Liouville 
description of a fluid assuming that microscopic fluctuations are linear and the system has 
no nonequilibrium fields. Green-Kubo formulae also assume the validity of the Onsager 
hypothesis, i.e., that spontaneous fluctuations in microscopic quantities decay according to 
hydrodynamic laws, and that hydrodynamic quantities are meaningful. This requires that 
time scales are long compared to the collision time and that the system size is large compared 
to the mean free path. 

To compute the shear viscosity, we start with time series data for the positions {x^ yi) 
and velocities (v X) i,v yi i) of N particles, as well as the shear stress 

P xy (t) = E mv x>i v y>i - E E ■ r ' j!J ' jl " inj) . (1) 

i=l i j>i ^"ij 

The first term of Eq. (0) is a kinetic part, which depends only on particle velocities, and the 
second term is a potential part, which depends on the pair potential. Here m is the particle 
mass and = (xij, yij) is the distance between particles % and j. We can then compute the 
shear stress autocorrelation function (SACF) 

C shear (t) = (P^(t)P^(O)). (2) 

Finally, we find r? by integrating the SACF using the Green-Kubo relation 

1 r°° 

V= Ak^fJ C ^(t)dt, (3) 

for a 2D liquid of area A and temperature T. Equation (j2J) yields the hydrodynamic param- 
eter t] based on fluctuating microscopic parameters entering into the shear stress P xy {t). 

We use normalized units in this letter. The length and time are normalized by a and 

n 

uj~J, respectively, where uj pd = (Q 2 /2ire ma 3 ) 1 / 2 [9]. The normalized temperature is r -1 , 
where T = Q 2 / 47ie akT is the Coulomb coupling parameter, so that a high temperature 
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corresponds to a small T. The 2D viscosity rj is normalized by 770 = nmu p da 2 , and the 
kinematic viscosity v = rj/nm is normalized by w p d,a 2 . 

We performed an MD simulation to calculate rj. The equations of motion for iV particles 
were integrated using periodic boundary conditions. A thermostat was applied to achieve 
a constant T. We recorded particle positions and velocities, and we used Eqs. (HJ)-© to 
calculate n. 

Our —on .node, resent the expe rim ental 8y8tem in Ref. Q. In both of ttan, 
particles in a monolayer interact with a Yukawa potential. The values of k and V were 
similar; all our simulations were performed for k = 0.56 while the experiment had k = 0.53. 
There are, however, significant differences. The simulation is for equilibrium conditions, 



while the experiment of Ref. 



12j, like most experiments to measure 77, used an externally- 



applied shear and therefore resulted in a measurement of 77 under nonequilibrium conditions. 
The simulation had periodic boundary conditions, unlike the experiment, and the equation 
of motion when a thermostat is used does not explicitly model the frictional damping of 
particle motion due to gas in the experiment. 

We now review the details and tests of our simulation method. We used a velocity Verlet 
integrator ^| with a time step 0.02 < At < 0.05 uj~J. We verified that At was small enough 
by performing a test, with no thermostat, where we required a fluctuation of total energy 
< 3% over an interval of 750 coQ 1 . We truncated the Yukawa potential at r cut = 22a, with 
a switching function to give a smooth cutoff between 20a < r < 22a. We verified that the 
potential energy of the entire system was almost independent of r cut , for r cut > 12a. We 
used N = 1024 particles, corresponding to a rectangle 56.99a x 49.08a. The size of this 
simulation box limits the maximum meaningful time for correlation functions to 46 u;"^ 1 , 
computed as the time for a compressional sound wave to transit the box. Later we will find 
that except for T > 124, which is near freezing, the SACF decays to zero in less than 46 
ujpdi indicating that our simulation box was chosen sufficient small. Ewald summation was 
not used because the simulation box was wider than A 73 by a factor of 27. The ratio of the 
two sides of the box were chosen to allow a perfect triangular lattice to form at high T, i.e., 
at low temperature. 

After completing these tests, we added the Nose-Hoover thermostat (2^ to the equation 
of motion. We tested the thermostat with different values of the thermal relaxation time, 
and we chose a value of 1.0 uj~ d , which resulted in a canonical distribution within a time 
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4000 Lo Z j . To verify that a thermal equilibrium was attained, we performed the customary 
test 21] of temperature fluctuations, characterized by their variance and skewness. We also 
verified that energy was equally partitioned among collective modes. 

Our results were prepared in four steps. First, an initial configuration of random particle 
positions and velocities was chosen. After T reached the desired level and equilibrium was 
attained, we began recording data for a duration of 4 x lO 4 ^/^ 1 . Second, the SACF was 
calculated using Eqs. (JIJ and (j2J) for the entire duration. To verify the validity of this result, 
we tested ten ensemble averages, each for a time series of a different duration, and we found 
that results for the SACF were independent of the duration, for the duration we recorded. 
Third, we integrated the SACF over t, yielding a value for the shear viscosity rj. To verify 
that 7] does not depend on N, we calculated 77 using N = 4096 for T = 17 and 124, and 
we found that rj was the same, within error bars, as for N = 1024. Fourth, we averaged 
the results for 3 to 6 different initial configurations, yielding our chief results the SACF, as 
shown in Fig. ^ and 77, as shown in Fig. Ufa). 

The SACF decays rapidly with time for T < 124. This decay is almost exponential with 
t for large T, i.e., \nC s h e ar{t) oc —t. It decays even faster, \nC s h e ar{t) oc —t 2 , for small T. 



These results are s 
previous results 



lown for T = 17 and 89 in Fig. H Our results are contrary to some 
for 2D systems, where a C s h ea r{t) oc t 1 dependence was found in 
the tail of the stress autocorrelation function. We will discuss this at the end of this letter. 

We find that the shear viscosity calculated using the Green-Kubo relation is finite in 2D 
liquids with a Yukawa potential. This is true for a wide range of T, except possibly near 
the freezing region. This is indicated by the exponential (or faster than exponential) decay 
of the SACF with time in Fig. so that the time integral of the SACF converges when the 
Green-Kubo relation is used to calculate 77. However, when the system is near freezing, i.e., 
r > 124, the decay is slower than exponential; thus, our result for 77 in this regime is less 
reliable. 

Our chief result is the variation of 77 with T, Fig. Ufa). At high temperature, 77 decreases 
with T. In this regime, the system behaves as a kinetic gas, corresponding to a disordered 
state, as seen from the orbits in Fig. Efb). When V is larger than 17, on the other hand, 
7] increases with T; it exhibits an exponential dependence on V for V < 124 and a much 
steeper increase for V > 124. Near freezing, T > 124, the system has a highly ordered 
structure, as seen from the orbits in Fig. E[c). The minimum viscosity, which is at T = 17, 
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is O.I4770. Using experimental parameters j3| a = 0.6 mm and u p d = 40 s _1 , our minimum 
corresponds to a kinematic viscosity v — 2 mm 2 s _1 in physical units. This is about two 
times larger than the kinematic viscosity of liquid water at STP conditions. This result, 



as noted previously 



24j, is true even though rj is itself an extremely small number for a 



dusty plasma. The reason is that the mass density is also very small, so that the kinematic 
viscosity, which is the ratio of these two small quantities, happens to be comparable to that 
of a denser substance like water. 

We compare our results with the experiment of Ref. in Fig. El As in our simulation, 
the experimental 77 varies with T, and it has a minimum. The minimum value rj = 0. 13?7o in 
the experiment matches our result of 0.1477 . 

Aside from this agreement in the magnitude of 77, however, there is a difference in the 
value of T where the minimum of 7/ occurs. The experimental result exhibits a much broader 
minimum, and the minimum occurs at a much higher T, as seen in Fig. El We suggest two 
reasons for these differences. Both of these reasons arise from inhomogeneity and anisotropy 
in the experiment that are lacking in the simulation. First, the experiment was nonequi- 
librium, with an applied shear that had a specific scale length and that was in a specific 
direction. In contrast, the simulation was in equilibrium, with the shear corresponding to 
thermal motions that had a wide range of length scales including very short length scales, 
and the direction of the shear fluctuated isotropically. Second, the experiment had a nonuni- 
form temperature; therefore it had nonuniform values of T and 7/ whereas the simulation 
was uniform. The values reported for 77 and T in experiment ^] were computed as spatial 
averages over a region that had a nonuniform temperature; this probably had its most sig- 
nificant effect on the value of T. Thus, it is not surprising that the T for the minimum in 
the experiment does not match that of the simulation. 

To discover the effect of the dimensionality of a system, we compare the kinematic vis- 
cosity of 2D and 3D liquids. Saigo and Hamaguchi |l5j| performed an equilibrium simulation 
similar to ours, with a Yukawa potential, and their results for k = 0.5 are plotted in Fig. 01 
In both cases the viscosity has a minimum at T « 20, but the magnitudes are not the same. 
In 2D, the kinematic viscosity is mostly larger for the same value of T. 

The minimum of 77 with temperature is a distinctive feature not found in most simple 
liquids. In water, for example, viscosity decreases monotonically with temperature. Systems 
such as strongly-coupled plasmas with a long-range repulsive potential, however, tend to 
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have a minimum. This minimum has been found in a 2D dusty plasma experiment 12 1, 



14. 



15 



simulations of liquids with Yukawa potential in 2D (the present work) and 3D 
and a simulation of a one-component plasma (OCP) [lsI ]. 

The minimum arises from the temperature dependence of the kinetic and potential con- 
tributions to momentum transport. This is seen in Fig. [2(d), where the kinetic part of 77 
decreases with T, while the potential part increases with T. (There is also a third contribu- 
tion, called the cross term [18j . but we found it is insignificant for our conditions.) Neither 
simple liquids nor dilute gases have a minimum because they are dominated by the potential 
and kinetic contributions, respectively. 

Finally, we discuss a controversy for shear viscosity in 2D liquids. Previous theoretical 
and simulation efforts, with a non- Yukawa potential, yielded conflicting results for the decay 
of the SACF. This is important because using Green-Kubo relations to compute transport 
coefficients requires a decay rapid enough for the integral to converge. Previous efforts 
using hydrodynamic mode-coupling theory (3] and an MD simulation Q] predicted a 
dependence for the tail in the SACF; this slow decay could lead to a divergent result for 
the viscosity. However, other MD simulations [u| yielded a much more rapid decay, which 
allows the Green-Kubo integral to converge. Our result for a Yukawa potential is consistent 
with the latter, not the former result. This is true for T < 124; for T > 124, our data for the 
decay was not conclusive, so that further study with a larger simulation box and more initial 
conditions is needed to resolve this controversy in that range, which is very near freezing. 

We thank V. Nosenko and F. Skiff for helpful discussions. This work was supported by 
NASA and DOE. 
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FIG. 1: Shear stress autocorrelation function C s / iear (t), computed using Eq. (J2J). It is significant 
that this function decays exponentially, or even faster, so that integrating it over time t, Eq. (J2J), 
yields a meaningful value of rj. 



FIG. 2: (a) Simulation results show that the 2D shear viscosity r\ varies with temperature T -1 , and 
it has a minimum at T = 17. Data are shown for k = 0.56. Error bars represent the uncertainty, 
(b) Trajectories of particles during a time interval of 1.0 uj~J for T = 1 indicate a state of complete 
disorder, (c) Trajectories of particles during a time interval of 1.0 uj~J for T = 124 indicate a highly 
ordered structure, (d) The shear viscosity r] is primarily the sum of two contributions kinetic and 
potential. 



FIG. 3: Comparison of our 2D simulation with a 2D experiment and a 3D simulation. Nosenko 
and Goree's experiment Q| used a 2D dusty plasma at k = 0.53 with an externally-applied velocity 
shear. Saigo and Hamaguchi's simulation [15| used a 3D liquid with a Yukawa potential at k = 0.50, 
in the absence of externally-applied shear. In all three cases, rj has a minimum. 
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